Pylops - seismic interpolation

Author: M.Ravasi

In this notebook I will implement and compare different ways to perform interpolation of seismic data onto a regular grid from data acquired by an irregularly sampled array of receivers.

The problem of seismic data interpolation is a very simple one, yet very ill-posed and hard to solve. Our forward modelling operator is a simple restriction or sampling operator which is applied along the spatial direction(s).

$$\mathbf{y} = \mathbf{R} \mathbf{x}$$

where $\mathbf{y} = [\mathbf{y}_{R1}^T, \mathbf{y}_{R2}^T,..., \mathbf{y}_{RN^T}]^T$ and each vector $\mathbf{y}_{Ri}$ contains all time samples recorded in the seismic data at the specific receiver $R_i$. Similarly, $\mathbf{x} = [\mathbf{x}_{r1}^T, \mathbf{x}_{r2}^T,..., \mathbf{x}_{rM}^T]$, contains all traces at the regularly and finely sampled idea array of receivers $r_i$.

The forward model simply selects a subset of sources/receivers from the input vector which is assumed to have a fine and regularly sampled array of receivers (which is what we wish to recover via inversion), while its adjoint places those values at the correct locations in the input vector.

In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

from scipy.signal import convolve

from pylops.utils                      import dottest
from pylops.utils.wavelets             import *
from pylops.utils.seismicevents        import *
from pylops.basicoperators             import *
from pylops.signalprocessing           import *
from pylops.waveeqprocessing           import *
from pylops.optimization.leastsquares  import *
from pylops.optimization.sparsity  import FISTA

2D

Linear events

In [2]:
###### Input parameters
par = {'ox':0,'dx':2,    'nx':140,
       'oy':0,'dy':2,    'ny':101,
       'ot':0,'dt':0.004,'nt':200,
       'f0': 20, 'nfmax': 210}

v       = 1500
t0_m    = [0.2, 0.4, 0.5]
theta_m = [0, 30, -60]
phi_m   = [0]
amp_m   = [1., -2, 0.5]

# Create axis
taxis, taxis2, xaxis, yaxis = makeaxis(par)

# Create wavelet
wav = ricker(taxis[:41], f0=par['f0'])[0]

# Generate model
m, x =  linear2d(xaxis, taxis, v, t0_m, theta_m, amp_m, wav)

# subsampling locations
perc_subsampling=0.4
Nsub=int(np.round(par['nx']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['nx']))[:Nsub])

# Create restriction operator
Rop = Restriction(par['nx']*par['nt'], iava, dims=(par['nx'], par['nt']), dir=0, dtype='float64')
dottest(Rop,Nsub*par['nt'],par['nx']*par['nt'], verb=True)

y     = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, par['nt'])
xadj = xadj.reshape( par['nx'], par['nt'])

# Apply mask
ymask = Rop.mask(x.flatten())

Dot test passed, v^T(Opu)=-174.786015 - u^T(Op^Tv)=-174.786015

Inversion in time-space domain with spatial smoothing

In [3]:
xinvf, _, _ = SeismicInterpolation(y, par['nx'], iava, kind='spatial', dottest=True,
                                   **dict(epsRs=[np.sqrt(0.1)], returninfo=False,
                                          damp=np.sqrt(1e-4), iter_lim=200 ,show=0))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=2.532455 - u^T(Op^Tv)=2.532455

Inversion in f-k domain with sparsity-promoting solver

In [4]:
# FFT of full data
nfft=2**10
FFTop = FFT2D(dims=[par['nx'], par['nt']], nffts=[nfft, nfft], sampling=[par['dx'], par['dt']])
p = FFTop*x.flatten()
p = np.reshape(p, (nfft, nfft))

# Inverse
xinvf, pinvf, cost = SeismicInterpolation(y, par['nx'], iava, kind='fk', nffts=(nfft, nfft), 
                                          sampling=(par['dx'], par['dt']), dottest=False,
                                          **dict(niter=50, eps=1e-2, returninfo=True))

fig, axs = plt.subplots(1, 6, figsize=(20, 6))
axs[0].imshow(np.fft.fftshift(np.abs(p[:, :nfft//2-1]),axes=0).T, 
              extent=(np.fft.fftshift(FFTop.f1)[0], np.fft.fftshift(FFTop.f1)[-1], 
                      FFTop.f2[nfft//2-1], FFTop.f2[0]))
axs[0].set_title('Original FFT')
axs[0].axis('tight')
axs[0].set_ylim(50, 0)
axs[0].set_xlim(-0.05, 0.05)
axs[1].imshow(np.fft.fftshift(np.abs(pinvf[:, :nfft//2-1]),axes=0).T,
              extent=(np.fft.fftshift(FFTop.f1)[0], np.fft.fftshift(FFTop.f1)[-1], 
                      FFTop.f2[nfft//2-1], FFTop.f2[0]))
axs[1].set_title('Reconstructed FFT')
axs[1].axis('tight')
axs[1].set_ylim(50, 0)
axs[1].set_xlim(-0.05, 0.05)
axs[2].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Data')
axs[2].axis('tight');
axs[3].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Masked Data')
axs[3].axis('tight');
axs[4].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Reconstructed data (sparse inv.)')
axs[4].axis('tight');
axs[5].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[5].set_title('Error')
axs[5].axis('tight');

Inversion in linear Radon domain with sparsity-promoting solver

First we just setup and apply Radon to the full data

In [5]:
npx=200
pxmax = 1e-3
px = np.linspace(-pxmax, pxmax, npx)
Radop = Radon2D(taxis, xaxis, px, engine='numba')
dottest(Radop, par['nx']*par['nt'], npx*par['nt'], verb=True)

Rop1 = Rop*Radop

# Radon on original data
xradadj = Radop.H*x.flatten()
xradadj = xradadj.reshape(npx, par['nt'])

yradadj = Radop*xradadj.flatten()
yradadj = yradadj.reshape(par['nx'], par['nt'])

# Adjoint
xadj = Rop1.H*y.flatten()
xadj = xadj.reshape(npx, par['nt'])

fig, axs = plt.subplots(1, 4, figsize=(20, 6))
axs[0].imshow(x.T, cmap='gray')
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj.T, cmap='gray')
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj.T, cmap='gray')
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj.T, cmap='gray')
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');

Dot test passed, v^T(Opu)=744.126259 - u^T(Op^Tv)=744.126259

Then we apply sparse inversion

In [6]:
xinvf, pinvf, cost = SeismicInterpolation(y, par['nx'], iava, kind='radon-linear', spataxis=xaxis, 
                                          taxis=taxis, paxis=px, centeredh=True, dottest=False,
                                          **dict(niter=20, eps=5e-1, returninfo=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(pinvf.T, cmap='gray')
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

Parabolic events

In [7]:
###### Input parameters
par = {'ox':0, 'dx':2,    'nx':140,
       'oy':0, 'dy':2,    'ny':100,
       'ot':0, 'dt':0.004,'nt':200,
       'f0': 20, 'nfmax': 210}

v       = 1500
t0_m    = [0.2, 0.4, 0.5]
px_m = [0, 0, 0]
pxx_m = [1e-5, 5e-6, 1e-6]
amp_m   = [1., -2, 0.5]

# Create axis
taxis, taxis2, xaxis, yaxis = makeaxis(par)

# Create wavelet
wav = ricker(taxis[:41], f0=par['f0'])[0]

# Generate model
_, x =  parabolic2d(xaxis, taxis, t0_m, px_m, pxx_m,amp_m,wav)

# subsampling locations
perc_subsampling=0.4
Nsub=int(np.round(par['nx']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['nx']))[:Nsub])

# Create restriction operator
Rop = Restriction(par['nx']*par['nt'], iava, dims=(par['nx'], par['nt']), dir=0, dtype='float64')
dottest(Rop,Nsub*par['nt'],par['nx']*par['nt'], verb=True)

y     = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, par['nt'])
xadj = xadj.reshape( par['nx'], par['nt'])

# Apply mask
ymask = Rop.mask(x.flatten())

Dot test passed, v^T(Opu)=10.398565 - u^T(Op^Tv)=10.398565

Inversion in time-space domain with spatial smoothing

In [8]:
xinvf, _, _ = SeismicInterpolation(y, par['nx'], iava, kind='spatial', dottest=True,
                                   **dict(epsRs=[np.sqrt(0.1)], returninfo=False,
                                          damp=np.sqrt(1e-4), iter_lim=200 ,show=0))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=123.862997 - u^T(Op^Tv)=123.862997

Inversion in f-k domain with sparsity-promoting solver

In [9]:
# FFT of full data
nfft=2**9
FFTop = FFT2D(dims=[par['nx'], par['nt']], nffts=[nfft, nfft], sampling=[par['dx'], par['dt']])
p = FFTop*x.flatten()
p = np.reshape(p, (nfft, nfft))

# Inverse
xinvf, pinvf, cost = SeismicInterpolation(y, par['nx'], iava, kind='fk', nffts=(nfft, nfft), 
                                          sampling=(par['dx'], par['dt']), dottest=True,
                                          **dict(niter=50, eps=1e-2, returninfo=True))

fig, axs = plt.subplots(1, 6, figsize=(20, 6))
axs[0].imshow(np.fft.fftshift(np.abs(p[:, :nfft//2-1]),axes=0).T, 
              extent=(np.fft.fftshift(FFTop.f1)[0], np.fft.fftshift(FFTop.f1)[-1], 
                      FFTop.f2[nfft//2-1], FFTop.f2[0]))
axs[0].set_title('Original FFT')
axs[0].axis('tight')
axs[0].set_ylim(50, 0)
axs[0].set_xlim(-0.05, 0.05)
axs[1].imshow(np.fft.fftshift(np.abs(pinvf[:, :nfft//2-1]),axes=0).T,
              extent=(np.fft.fftshift(FFTop.f1)[0], np.fft.fftshift(FFTop.f1)[-1], 
                      FFTop.f2[nfft//2-1], FFTop.f2[0]))
axs[1].set_title('Reconstructed FFT')
axs[1].axis('tight')
axs[1].set_ylim(50, 0)
axs[1].set_xlim(-0.05, 0.05)
axs[2].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Data')
axs[2].axis('tight');
axs[3].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Masked Data')
axs[3].axis('tight');
axs[4].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Reconstructed data (sparse inv.)')
axs[4].axis('tight');
axs[5].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[5].set_title('Error')
axs[5].axis('tight');

plt.figure(figsize=(7, 2))
plt.plot(cost, 'k', lw=2);

Dot test passed, v^T(Opu)=12.807560 - u^T(Op^Tv)=12.807560

Inversion in parabolic Radon domain with sparsity-promoting solver

First we just setup and apply Radon to the full data

In [10]:
npx=200
pxmax = 5e-5
px = np.linspace(-pxmax, pxmax, npx)
Radop = Radon2D(taxis, xaxis, px, centeredh=False, kind='parabolic', engine='numba')
dottest(Radop, par['nx']*par['nt'], npx*par['nt'], verb=True)

Rop1 = Rop*Radop

# Radon on original data
xradadj = Radop.H*x.flatten()
xradadj = xradadj.reshape(npx, par['nt'])

yradadj = Radop*xradadj.flatten()
yradadj = yradadj.reshape(par['nx'], par['nt'])

# Adjoint
xadj = Rop1.H*y.flatten()
xadj = xadj.reshape(npx, par['nt'])

fig, axs = plt.subplots(1, 4, figsize=(20, 6))
axs[0].imshow(x.T, cmap='gray')
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj.T, cmap='gray')
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj.T, cmap='gray')
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj.T, cmap='gray')
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');

Dot test passed, v^T(Opu)=606.248786 - u^T(Op^Tv)=606.248786

Then we apply sparse inversion

In [11]:
xinvf, pinvf, _ = SeismicInterpolation(y, par['nx'], iava, kind='radon-parabolic', spataxis=xaxis, 
                                       taxis=taxis, paxis=px, centeredh=False, dottest=False,
                                       **dict(niter=20, eps=5e-1, returninfo=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(pinvf.T, cmap='gray')
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

Inversion in sliding linear Radon domain with sparsity-promoting solver

In [12]:
nwin = 26
nwins = 8
nover = 10
npx = 61
px = np.linspace(-5e-3, 5e-3, npx)
dimsd = x.shape
dims = (nwins*npx, dimsd[1])

# sliding window radon with overlap
Op = Radon2D(taxis, np.linspace(-par['dx']*nwin//2, par['dx']*nwin//2, nwin), px, centeredh=True,
             kind='linear', engine='numba')
Slidop = Sliding2D(Op, dims, dimsd, nwin, nover, tapertype='cosine')
dottest(Slidop, np.prod(dimsd), np.prod(dims))

Rop1 = Rop*Slidop

# adjoint
padj = Rop1.H*y.flatten()
padj = padj.reshape(npx*nwins, par['nt'])

# inverse
xinvf, pinvf, _ = SeismicInterpolation(y, par['nx'], iava, kind='sliding', spataxis=xaxis, 
                                    taxis=taxis, paxis=px, nwins=nwins, nwin=nwin, 
                                    nover=nover, dottest=False,
                                    **dict(niter=40, eps=7e-1, returninfo=True, show=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(pinvf.T, cmap='gray')
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x.T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask.T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x.T-xinvf.T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

FISTA optimization
-----------------------------------------------------------
The Operator Op has 11200 rows and 97600 cols
eps = 7.000000e-01	tol = 1.000000e-10	niter = 40
alpha = 1.668533e-03	thresh = 5.839865e-04
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1  -0.00000e+00   4.238e+02   4.512e+02   4.377e-01
     2  -0.00000e+00   3.407e+02   3.902e+02   3.615e-01
     3  -0.00000e+00   2.648e+02   3.373e+02   3.901e-01
     4  -0.00000e+00   2.015e+02   2.961e+02   3.928e-01
     5  -0.00000e+00   1.519e+02   2.665e+02   3.782e-01
     6  -0.00000e+00   1.147e+02   2.463e+02   3.532e-01
     7  -0.00000e+00   8.760e+01   2.330e+02   3.235e-01
     8  -0.00000e+00   6.818e+01   2.241e+02   2.927e-01
     9  -0.00000e+00   5.444e+01   2.180e+02   2.634e-01
    10  -0.00000e+00   4.479e+01   2.132e+02   2.373e-01
    11  -0.00000e+00   3.804e+01   2.091e+02   2.149e-01
    21   0.00000e+00   2.299e+01   1.725e+02   1.324e-01
    31   0.00000e+00   2.182e+01   1.551e+02   9.906e-02
    32   0.00000e+00   2.173e+01   1.542e+02   9.554e-02
    33   0.00000e+00   2.164e+01   1.534e+02   9.243e-02
    34   0.00000e+00   2.155e+01   1.526e+02   8.956e-02
    35   0.00000e+00   2.147e+01   1.520e+02   8.622e-02
    36   0.00000e+00   2.139e+01   1.514e+02   8.318e-02
    37   0.00000e+00   2.132e+01   1.509e+02   8.085e-02
    38   0.00000e+00   2.125e+01   1.503e+02   7.869e-02
    39  -0.00000e+00   2.119e+01   1.498e+02   7.660e-02
    40  -0.00000e+00   2.113e+01   1.494e+02   7.462e-02

FISTA finished

Realistic data

In [13]:
inputfile = '../../pylops/testdata/marchenko/input.npz'
inputdata = np.load(inputfile)

x = inputdata['R'][50, :, :-100]
x = x/np.abs(x).max()

taxis, xaxis = inputdata['t'][:-100], inputdata['r'][0]
par = {}
par['nx'], par['nt'] = x.shape
par['dx'] = inputdata['r'][0, 1] - inputdata['r'][0, 0]
par['dt'] = inputdata['t'][1] - inputdata['t'][0]

# add wavelet
wav = 0.3*inputdata['wav']
wav_c = np.argmax(wav)
x = np.apply_along_axis(convolve, 1, x, wav, mode='full')
x = x[:, wav_c:][:, :par['nt']]

# gain
gain = np.tile((taxis**2)[:, np.newaxis], (1, par['nx'])).T
x = x*gain

# subsampling locations
perc_subsampling=0.6
Nsub=int(np.round(par['nx']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['nx']))[:Nsub])

# Create restriction operator
Rop = Restriction(par['nx']*par['nt'], iava, dims=(par['nx'], par['nt']), dir=0, dtype='float64')
dottest(Rop,Nsub*par['nt'],par['nx']*par['nt'], verb=True)

y     = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, par['nt'])
xadj = xadj.reshape( par['nx'], par['nt'])

# Apply mask
ymask = Rop.mask(x.flatten())

Dot test passed, v^T(Opu)=-45.086423 - u^T(Op^Tv)=-45.086423

Inversion in time-space domain with spatial smoothing

In [14]:
xinv, _, _ = SeismicInterpolation(y, par['nx'], iava, kind='spatial', dottest=True,
                               **dict(epsRs=[np.sqrt(0.1)], returninfo=False,
                                      damp=np.sqrt(1e-4), iter_lim=200, show=0))

fig, axs = plt.subplots(1, 5, figsize=(20, 9))
axs[0].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x.T-xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=66.988570 - u^T(Op^Tv)=66.988570

Inversion in f-k domain with sparsity-promoting solver

In [15]:
# FFT of full data
nfft=2**10
FFTop = FFT2D(dims=[par['nx'], par['nt']], nffts=[nfft, nfft], sampling=[par['dx'], par['dt']])
p = FFTop*x.flatten()
p = np.reshape(p, (nfft, nfft))

pmask = FFTop*x.flatten()
pmask = np.reshape(p, (nfft, nfft))


# Inverse
xinv, pinv, _ = SeismicInterpolation(y, par['nx'], iava, kind='fk', nffts=(nfft, nfft), 
                                     sampling=(par['dx'], par['dt']), dottest=False,
                                     **dict(niter=100, eps=1e-6, returninfo=True, show=False))

fig, axs = plt.subplots(1, 5, figsize=(20, 9))
axs[0].imshow(np.fft.fftshift(np.abs(p[:, :nfft//2-1]),axes=0).T, 
              extent=(np.fft.fftshift(FFTop.f1)[0], np.fft.fftshift(FFTop.f1)[-1], 
                      FFTop.f2[nfft//2-1], FFTop.f2[0]))
axs[0].set_title('Original FFT')
axs[0].axis('tight')
axs[0].set_ylim(70, 0)
axs[1].imshow(np.fft.fftshift(np.abs(pinv[:, :nfft//2-1]),axes=0).T,
              extent=(np.fft.fftshift(FFTop.f1)[0], np.fft.fftshift(FFTop.f1)[-1], 
                      FFTop.f2[nfft//2-1], FFTop.f2[0]))
axs[1].set_title('Reconstructed FFT')
axs[1].axis('tight')
axs[1].set_ylim(70, 0)
axs[2].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[2].set_title('Data')
axs[2].axis('tight');
axs[3].imshow(xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x.T-xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[4].set_title('Error')
axs[4].axis('tight');

Inversion in parabolic Radon domain with sparsity-promoting solver

In [16]:
npx=201
pxmax = 4e-5
px = np.linspace(-pxmax/3, pxmax, npx)
Radop = Radon2D(taxis, xaxis, px, centeredh=True, kind='parabolic', engine='numba')
dottest(Radop, par['nx']*par['nt'], npx*par['nt'], verb=True)

Rop1 = Rop*Radop

# Radon on original data
xradadj = Radop.H*x.flatten()
xradadj = xradadj.reshape(npx, par['nt'])

yradadj = Radop*xradadj.flatten()
yradadj = yradadj.reshape(par['nx'], par['nt'])

# Adjoint
xadj = Rop1.H*y.flatten()
xadj = xadj.reshape(npx, par['nt'])

fig, axs = plt.subplots(1, 4, figsize=(20, 9))
axs[0].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj.T, cmap='gray', vmin=-1, vmax=1)
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj.T, cmap='gray', vmin=-50, vmax=50)
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj.T, cmap='gray', vmin=-1, vmax=1)
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');

# Sparse inversion
xinv, pinv, _ = SeismicInterpolation(y, par['nx'], iava, kind='radon-parabolic', spataxis=xaxis, 
                                     taxis=taxis, paxis=px, centeredh=True, dottest=False,
                                     **dict(niter=50, eps=4e-2, alpha=1e-4, returninfo=True, show=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 9))
axs[0].imshow(pinv.T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x.T-xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=-1230.283443 - u^T(Op^Tv)=-1230.283443
FISTA optimization
-----------------------------------------------------------
The Operator Op has 42700 rows and 140700 cols
eps = 4.000000e-02	tol = 1.000000e-10	niter = 50
alpha = 1.000000e-04	thresh = 2.000000e-06
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   0.00000e+00   8.863e-01   8.998e-01   2.238e-03
     2   0.00000e+00   8.362e-01   8.624e-01   2.119e-03
     3   0.00000e+00   7.785e-01   8.198e-01   2.577e-03
     4   0.00000e+00   7.164e-01   7.748e-01   2.949e-03
     5   0.00000e+00   6.529e-01   7.296e-01   3.242e-03
     6   0.00000e+00   5.902e-01   6.859e-01   3.463e-03
     7   0.00000e+00   5.299e-01   6.450e-01   3.621e-03
     8   0.00000e+00   4.733e-01   6.076e-01   3.723e-03
     9   0.00000e+00   4.209e-01   5.739e-01   3.779e-03
    10   0.00000e+00   3.731e-01   5.440e-01   3.796e-03
    11   0.00000e+00   3.298e-01   5.177e-01   3.782e-03
    21   0.00000e+00   9.402e-02   3.826e-01   2.949e-03
    31   0.00000e+00   3.734e-02   3.263e-01   2.163e-03
    41   0.00000e+00   2.714e-02   2.850e-01   1.720e-03
    42   0.00000e+00   2.681e-02   2.816e-01   1.686e-03
    43   0.00000e+00   2.654e-02   2.784e-01   1.654e-03
    44   0.00000e+00   2.631e-02   2.753e-01   1.623e-03
    45   0.00000e+00   2.611e-02   2.723e-01   1.592e-03
    46   0.00000e+00   2.595e-02   2.695e-01   1.562e-03
    47   0.00000e+00   2.580e-02   2.668e-01   1.534e-03
    48   0.00000e+00   2.568e-02   2.643e-01   1.507e-03
    49   0.00000e+00   2.557e-02   2.619e-01   1.480e-03
    50   0.00000e+00   2.547e-02   2.596e-01   1.450e-03

FISTA finished

Inversion in hyperbolic Radon domain with sparsity-promoting solver

In [17]:
npx=201
pxmax = 5e-5
px = np.linspace(0, pxmax, npx)
Radop = Radon2D(taxis, xaxis, px, centeredh=True, kind='hyperbolic', engine='numba')
dottest(Radop, par['nx']*par['nt'], npx*par['nt'], verb=True)

Rop1 = Rop*Radop

# Radon on original data
xradadj = Radop.H*x.flatten()
xradadj = xradadj.reshape(npx, par['nt'])

yradadj = Radop*xradadj.flatten()
yradadj = yradadj.reshape(par['nx'], par['nt'])

# Adjoint
xadj = Rop1.H*y.flatten()
xadj = xadj.reshape(npx, par['nt'])

fig, axs = plt.subplots(1, 4, figsize=(20, 9))
axs[0].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj.T, cmap='gray', vmin=-1, vmax=1)
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj.T, cmap='gray', vmin=-50, vmax=50)
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj.T, cmap='gray', vmin=-1, vmax=1)
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');

# Sparse inversion
xinv, pinv, _ = SeismicInterpolation(y, par['nx'], iava, kind='radon-hyperbolic', spataxis=xaxis, 
                                     taxis=taxis, paxis=px, centeredh=True, dottest=False,
                                     **dict(niter=50, eps=4e-2, returninfo=True, show=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 9))
axs[0].imshow(pinv.T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x.T-xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=3102.490485 - u^T(Op^Tv)=3102.490485
FISTA optimization
-----------------------------------------------------------
The Operator Op has 42700 rows and 140700 cols
eps = 4.000000e-02	tol = 1.000000e-10	niter = 50
alpha = 8.541950e-05	thresh = 1.708390e-06
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   0.00000e+00   7.683e-01   7.917e-01   3.947e-03
     2   0.00000e+00   6.574e-01   6.983e-01   3.099e-03
     3   0.00000e+00   5.595e-01   6.180e-01   3.281e-03
     4   0.00000e+00   4.765e-01   5.514e-01   3.331e-03
     5   0.00000e+00   4.060e-01   4.959e-01   3.333e-03
     6   0.00000e+00   3.456e-01   4.488e-01   3.319e-03
     7   0.00000e+00   2.936e-01   4.085e-01   3.293e-03
     8   0.00000e+00   2.489e-01   3.741e-01   3.250e-03
     9   0.00000e+00   2.108e-01   3.450e-01   3.184e-03
    10   0.00000e+00   1.787e-01   3.207e-01   3.097e-03
    11   0.00000e+00   1.518e-01   3.006e-01   2.991e-03
    21   0.00000e+00   5.406e-02   2.278e-01   1.569e-03
    31   0.00000e+00   4.747e-02   2.132e-01   8.497e-04
    41   0.00000e+00   4.607e-02   2.038e-01   6.847e-04
    42   0.00000e+00   4.595e-02   2.031e-01   6.734e-04
    43   0.00000e+00   4.582e-02   2.026e-01   6.623e-04
    44   0.00000e+00   4.570e-02   2.020e-01   6.516e-04
    45   0.00000e+00   4.557e-02   2.015e-01   6.415e-04
    46   0.00000e+00   4.543e-02   2.011e-01   6.311e-04
    47   0.00000e+00   4.529e-02   2.006e-01   6.212e-04
    48   0.00000e+00   4.515e-02   2.002e-01   6.118e-04
    49   0.00000e+00   4.501e-02   1.999e-01   6.028e-04
    50   0.00000e+00   4.487e-02   1.995e-01   5.940e-04

FISTA finished

Inversion in sliding linear Radon domain with sparsity-promoting solver

In [18]:
dx = par['dx']
nwin=20
nwins=9
nover=10
npx=101
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)
dimsd = x.shape
dims = (nwins*npx, dimsd[1])

# sliding window radon with overlap
Op = Radon2D(taxis, np.linspace(-dx*nwin//2, dx*nwin//2, nwin), px, centeredh=True,
             kind='linear', engine='numba')
Slidop = Sliding2D(Op, dims, dimsd, nwin, nover, tapertype='cosine', design=False)
dottest(Slidop, np.prod(dimsd), np.prod(dims))

# radon on one patch
x1 = x[1*nwin:2*nwin]
y1 = Op.H*x1.flatten()
y1 = y1.reshape(npx, par['nt'])

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(5, 7))
axs[0].imshow(x1.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[0].axis('tight')
axs[0].set_title('Single patch')

axs[1].imshow(y1.T, cmap='gray', vmin=-0.5, vmax=0.5)
axs[1].axis('tight')
axs[1].set_title('Adjoint Radon of Single patch');
In [19]:
# adjoint
Rop1 = Rop*Slidop

padj_fromx = Slidop.H*x.flatten()
padj_fromx = padj_fromx.reshape(npx*nwins, par['nt'])

padj = Rop1.H*y.flatten()
padj = padj.reshape(npx*nwins, par['nt'])

# inverse
xinv, pinv, _ = SeismicInterpolation(y, par['nx'], iava, kind='sliding', spataxis=xaxis, 
                                     taxis=taxis, paxis=px, nwins=nwins, nwin=nwin, 
                                     nover=nover, dottest=True,
                                     **dict(niter=50, eps=5e-3, alpha=1.2e-3, 
                                            returninfo=True, show=True))

fig, axs = plt.subplots(1, 6, figsize=(25, 9))
axs[0].imshow(padj_fromx.T, cmap='gray', vmin=-0.5, vmax=0.5)
axs[0].set_title('Adjoint Radon on data')
axs[0].axis('tight')
axs[1].imshow(padj.T, cmap='gray', vmin=-0.5, vmax=0.5)
axs[1].set_title('Adjoint Radon on subsampled data')
axs[1].axis('tight')
axs[2].imshow(pinv.T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Inverse Radon on subsampled data')
axs[2].axis('tight')
axs[3].imshow(x.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[3].set_title('Data')
axs[3].axis('tight');
axs[4].imshow(xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[4].set_title('Reconstructed data (sparse inv.)')
axs[4].axis('tight')
axs[5].imshow(x.T-xinv.T, cmap='gray', vmin=-0.05, vmax=0.05)
axs[5].set_title('Error')
axs[5].axis('tight');

Dot test passed, v^T(Opu)=-986.142461 - u^T(Op^Tv)=-986.142461
FISTA optimization
-----------------------------------------------------------
The Operator Op has 42700 rows and 636300 cols
eps = 5.000000e-03	tol = 1.000000e-10	niter = 50
alpha = 1.200000e-03	thresh = 3.000000e-06
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   0.00000e+00   7.451e-01   7.709e-01   1.542e-02
     2   0.00000e+00   5.981e-01   6.458e-01   1.323e-02
     3   0.00000e+00   4.575e-01   5.289e-01   1.469e-02
     4   0.00000e+00   3.358e-01   4.314e-01   1.516e-02
     5   0.00000e+00   2.384e-01   3.569e-01   1.486e-02
     6   0.00000e+00   1.653e-01   3.046e-01   1.401e-02
     7   0.00000e+00   1.133e-01   2.705e-01   1.278e-02
     8   0.00000e+00   7.807e-02   2.500e-01   1.136e-02
     9   0.00000e+00   5.510e-02   2.385e-01   9.883e-03
    10   0.00000e+00   4.058e-02   2.323e-01   8.470e-03
    11   0.00000e+00   3.159e-02   2.286e-01   7.208e-03
    21   0.00000e+00   1.945e-02   1.931e-01   3.560e-03
    31   0.00000e+00   1.928e-02   1.711e-01   3.025e-03
    41   0.00000e+00   1.827e-02   1.604e-01   2.629e-03
    42   0.00000e+00   1.823e-02   1.596e-01   2.590e-03
    43   0.00000e+00   1.820e-02   1.589e-01   2.551e-03
    44   0.00000e+00   1.817e-02   1.581e-01   2.512e-03
    45   0.00000e+00   1.815e-02   1.574e-01   2.476e-03
    46   0.00000e+00   1.813e-02   1.568e-01   2.440e-03
    47   0.00000e+00   1.810e-02   1.561e-01   2.401e-03
    48   0.00000e+00   1.808e-02   1.556e-01   2.362e-03
    49   0.00000e+00   1.805e-02   1.550e-01   2.326e-03
    50   0.00000e+00   1.803e-02   1.545e-01   2.290e-03

FISTA finished

3D

Linear events

In [20]:
###### Input parameters
par = {'oy':0,'dy':2,    'ny':51,
       'ox':0,'dx':2,    'nx':30,
       'ot':0,'dt':0.004,'nt':100,
       'f0': 20, 'nfmax': 210}

v       = 1500
t0_m    = [0.1, 0.2, 0.3]
theta_m = [0, 30, -60]
phi_m   = [0, 50, 30]
amp_m   = [1., -2, 0.5]

# Create axis
taxis, taxis2, xaxis, yaxis = makeaxis(par)

# Create wavelet
wav = ricker(taxis[:41], f0=par['f0'])[0]

# Generate model
m, x =  linear3d(xaxis, yaxis, taxis, v, t0_m, theta_m, phi_m, amp_m, wav)

# subsampling locations
perc_subsampling=0.4
Nsub=int(np.round(par['ny']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['ny']))[:Nsub])

# Create restriction operator
Rop = Restriction(par['ny']*par['nx']*par['nt'], iava, 
                  dims=(par['ny'], par['nx'], par['nt']), dir=0, dtype='float64')
dottest(Rop,Nsub*par['nx']*par['nt'], par['ny']*par['nx']*par['nt'], verb=True)

y = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, par['nx'], par['nt'])
xadj = xadj.reshape( par['ny'], par['nx'], par['nt'])

# Apply mask
ymask = Rop.mask(x.flatten())

Dot test passed, v^T(Opu)=-97.322443 - u^T(Op^Tv)=-97.322443

Inversion in time-space domain with spatial smoothing

In [21]:
# Invert for interpolated signal
xinv, _, _ = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='spatial', dottest=True,
                                  **dict(epsRs=[np.sqrt(0.1)], returninfo=False,
                                         damp=np.sqrt(1e-4), iter_lim=100, show=0))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[par['ny']//2].T-xinv[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=463.443772 - u^T(Op^Tv)=463.443772

Inversion in f-k domain with sparsity-promoting solver

In [22]:
nfft=2**7
FFTop = FFTND(dims=[par['ny'], par['nx'], par['nt']], nffts=[nfft, nfft, nfft], 
              sampling=[par['dy'], par['dx'], par['dt']])
dottest(FFTop, nfft*nfft*nfft, par['ny']*par['nx']*par['nt'], complexflag=2, verb=True)
p = FFTop*x.flatten()
p = np.reshape(p, (nfft, nfft, nfft))

# Inverse
xinv, pinv, _ = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='fk', nffts=(nfft, nfft, nfft), 
                                     sampling=(par['dy'], par['dx'], par['dt']), dottest=False,
                                     **dict(niter=50, eps=1e-1, alpha=1e0, returninfo=True, show=True))

fig, axs = plt.subplots(1, 6, figsize=(20, 9))
axs[0].imshow(np.fft.fftshift(np.abs(p[:, 1, :nfft//2-1]),axes=0).T, 
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      FFTop.fs[2][nfft//2-1], FFTop.fs[2][0]), vmin=0, vmax=10)
axs[0].set_title('Original FFT')
axs[0].axis('tight')
axs[0].set_ylim(50, 0)
axs[0].set_xlim(-0.05, 0.05)
axs[1].imshow(np.fft.fftshift(np.abs(pinv[:, 1, :nfft//2-1]), axes=0).T,
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      FFTop.fs[2][nfft//2-1], FFTop.fs[2][0]), vmin=0, vmax=10)
axs[1].set_title('Reconstructed FFT')
axs[1].axis('tight')
axs[1].set_ylim(50, 0)
axs[1].set_xlim(-0.05, 0.05)
axs[2].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Data')
axs[2].axis('tight');
axs[3].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Masked data')
axs[3].axis('tight')
axs[4].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Reconstructed data (sparse inv.)')
axs[4].axis('tight');
axs[5].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[5].set_title('Error')
axs[5].axis('tight');

Dot test passed, v^T(Opu)=659.528795 - u^T(Op^Tv)=659.528795
FISTA optimization
-----------------------------------------------------------
The Operator Op has 60000 rows and 2097152 cols
eps = 1.000000e-01	tol = 1.000000e-10	niter = 50
alpha = 1.000000e+00	thresh = 5.000000e-02
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   0.00000e+00   1.184e+03   2.783e+03   9.860e+01
     2   0.00000e+00   4.913e+02   2.344e+03   2.831e+01
     3   0.00000e+00   2.412e+02   2.158e+03   1.900e+01
     4   0.00000e+00   1.572e+02   2.027e+03   1.368e+01
     5   0.00000e+00   1.269e+02   1.901e+03   1.184e+01
     6  -0.00000e+00   1.141e+02   1.780e+03   1.144e+01
     7  -0.00000e+00   1.071e+02   1.672e+03   1.124e+01
     8  -0.00000e+00   1.016e+02   1.580e+03   1.100e+01
     9  -0.00000e+00   9.659e+01   1.500e+03   1.080e+01
    10  -0.00000e+00   9.184e+01   1.432e+03   1.054e+01
    11  -0.00000e+00   8.766e+01   1.372e+03   1.027e+01
    21  -0.00000e+00   7.272e+01   1.030e+03   8.370e+00
    31  -0.00000e+00   6.873e+01   8.836e+02   7.068e+00
    41  -0.00000e+00   6.718e+01   8.124e+02   5.817e+00
    42  -0.00000e+00   6.696e+01   8.075e+02   5.700e+00
    43  -0.00000e+00   6.675e+01   8.029e+02   5.611e+00
    44  -0.00000e+00   6.658e+01   7.983e+02   5.556e+00
    45  -0.00000e+00   6.646e+01   7.935e+02   5.528e+00
    46  -0.00000e+00   6.642e+01   7.891e+02   5.479e+00
    47  -0.00000e+00   6.636e+01   7.851e+02   5.404e+00
    48  -0.00000e+00   6.630e+01   7.814e+02   5.332e+00
    49  -0.00000e+00   6.622e+01   7.779e+02   5.252e+00
    50  -0.00000e+00   6.617e+01   7.745e+02   5.184e+00

FISTA finished

Inversion in linear Radon domain with sparsity-promoting solver

In [23]:
npx=40
pxmax = 4e-3
px = np.linspace(-pxmax, pxmax, npx)
Radop = Radon3D(taxis, yaxis, xaxis, px, px, engine='numba', interp=False)
dottest(Radop, par['ny']*par['nx']*par['nt'], npx*npx*par['nt'], verb=True)

Rop1 = Rop*Radop

# Radon on original data
xradadj = Radop.H*x.flatten()
xradadj = xradadj.reshape(npx, npx, par['nt'])

yradadj = Radop*xradadj.flatten()
yradadj = yradadj.reshape(par['ny'], par['nx'], par['nt'])

# Adjoint
xadj = Rop1.H*y.flatten()
xadj = xadj.reshape(npx, npx, par['nt'])

fig, axs = plt.subplots(1, 4, figsize=(20, 6))
axs[0].imshow(x[:, par['nx']//2].T, cmap='gray')
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj[:, npx//2].T, cmap='gray')
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj[:, par['nx']//2].T, cmap='gray')
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj[:, npx//2].T, cmap='gray')
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');

Dot test passed, v^T(Opu)=-11394.425972 - u^T(Op^Tv)=-11394.425972
In [24]:
xinv, pinv, _ = SeismicInterpolation(y, (par['ny'], par['nx']), 
                                     iava, kind='radon-linear', spataxis=yaxis, spat1axis=xaxis, 
                                     taxis=taxis, paxis=px, p1axis=px, centeredh=True, dottest=False,
                                     **dict(niter=50, eps=5e1, alpha=1.36e-6, returninfo=True, show=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(pinv[:, npx//2].T, cmap='gray')
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

FISTA optimization
-----------------------------------------------------------
The Operator Op has 60000 rows and 160000 cols
eps = 5.000000e+01	tol = 1.000000e-10	niter = 50
alpha = 1.360000e-06	thresh = 3.400000e-05
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   0.00000e+00   5.815e+03   5.857e+03   1.202e-02
     2   0.00000e+00   5.694e+03   5.778e+03   1.173e-02
     3   0.00000e+00   5.548e+03   5.681e+03   1.469e-02
     4   0.00000e+00   5.379e+03   5.570e+03   1.737e-02
     5   0.00000e+00   5.193e+03   5.448e+03   1.980e-02
     6   0.00000e+00   4.993e+03   5.319e+03   2.199e-02
     7   0.00000e+00   4.784e+03   5.184e+03   2.396e-02
     8   0.00000e+00   4.569e+03   5.047e+03   2.572e-02
     9   0.00000e+00   4.351e+03   4.909e+03   2.726e-02
    10   0.00000e+00   4.133e+03   4.772e+03   2.861e-02
    11   0.00000e+00   3.916e+03   4.637e+03   2.977e-02
    21   0.00000e+00   2.114e+03   3.587e+03   3.353e-02
    31   0.00000e+00   1.129e+03   3.044e+03   3.013e-02
    41   0.00000e+00   6.714e+02   2.713e+03   2.589e-02
    42   0.00000e+00   6.429e+02   2.684e+03   2.551e-02
    43   0.00000e+00   6.165e+02   2.656e+03   2.514e-02
    44   0.00000e+00   5.922e+02   2.630e+03   2.477e-02
    45   0.00000e+00   5.698e+02   2.603e+03   2.440e-02
    46   0.00000e+00   5.491e+02   2.577e+03   2.406e-02
    47   0.00000e+00   5.300e+02   2.552e+03   2.372e-02
    48   0.00000e+00   5.123e+02   2.527e+03   2.339e-02
    49   0.00000e+00   4.959e+02   2.502e+03   2.309e-02
    50   0.00000e+00   4.807e+02   2.477e+03   2.281e-02

FISTA finished

Inversion in sliding linear Radon domain with sparsity-promoting solver

In [25]:
dy = par['dy']
dx = par['dx']
nwin=(24, 11)
nwins=(2, 3)
nover=(2, 3)
npy=41
npx=41
pxmax = 5e-3
py = px = np.linspace(-pxmax, pxmax, npx)
dimsd = (par['ny'], par['nx'], par['nt'])
dims = (nwins[0]*npy, nwins[1]*npx, dimsd[2])

# sliding window radon with overlap
Op = Radon3D(taxis, np.linspace(-dy*nwin[0]//2, dy*nwin[0]//2, nwin[0]), 
             np.linspace(-dx*nwin[1]//2, dx*nwin[1]//2, nwin[1]), 
             px, px, centeredh=True,
             kind='linear', engine='numba')
Slidop = Sliding3D(Op, dims, dimsd, nwin, nover, (npy, npx), tapertype='cosine', design=True)
dottest(Slidop, np.prod(dimsd), np.prod(dims), verb=True)

# adjoint
Rop1 = Rop*Slidop

padj_fromx = Slidop.H*x.flatten()
padj_fromx = padj_fromx.reshape(nwins[0], nwins[1], npy, npx, par['nt'])

padj = Rop1.H*y.flatten()
padj = padj.reshape(nwins[0], nwins[1], npy, npx, par['nt'])

WARNING: (2,3) windows required...
WARNING: model wins - start0:[ 0 41], end0:[41 82], start1:[ 0 41 82], end1:[ 41  82 123]
WARNING: data wins - start0:[ 0 22], end0:[24 46], start1:[ 0  8 16], end1:[11 19 27]

Dot test passed, v^T(Opu)=-7564.937654 - u^T(Op^Tv)=-7564.937654
In [26]:
xinvf, pinvf, _ = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='sliding', 
                                       spataxis=yaxis, spat1axis=xaxis, 
                                       taxis=taxis, paxis=py, p1axis=px, nwins=nwins, nwin=nwin, 
                                       nover=nover, dottest=False,
                                       **dict(niter=40, eps=5e0, alpha=1e-5, returninfo=True, show=True))

fig, axs = plt.subplots(2, 4, figsize=(25, 12))
axs[0][0].imshow(padj_fromx[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-100, vmax=100)
axs[0][0].set_title('Adjoint Radon on data')
axs[0][0].axis('tight')
axs[0][1].imshow(padj[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-100, vmax=100)
axs[0][1].set_title('Adjoint Radon on subsampled data')
axs[0][1].axis('tight')
axs[0][2].imshow(pinvf[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-0.1, vmax=0.1)
axs[0][2].set_title('Inverse Radon on subsampled data')
axs[0][2].axis('tight')
axs[1][0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1][0].set_title('Data')
axs[1][0].axis('tight');
axs[1][1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1][1].set_title('Masked data')
axs[1][1].axis('tight');
axs[1][2].imshow(xinvf[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1][2].set_title('Reconstructed data (sparse inv.)')
axs[1][2].axis('tight')
axs[1][3].imshow(x[:, par['nx']//2].T-xinvf[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1][3].set_title('Error')
axs[1][3].axis('tight');

FISTA optimization
-----------------------------------------------------------
The Operator Op has 60000 rows and 1008600 cols
eps = 5.000000e+00	tol = 1.000000e-10	niter = 40
alpha = 1.000000e-05	thresh = 2.500000e-05
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1  -0.00000e+00   5.531e+03   5.631e+03   6.202e-02
     2  -0.00000e+00   5.204e+03   5.387e+03   5.472e-02
     3  -0.00000e+00   4.857e+03   5.132e+03   6.379e-02
     4  -0.00000e+00   4.512e+03   4.883e+03   7.026e-02
     5  -0.00000e+00   4.183e+03   4.653e+03   7.456e-02
     6  -0.00000e+00   3.879e+03   4.447e+03   7.705e-02
     7  -0.00000e+00   3.605e+03   4.268e+03   7.809e-02
     8  -0.00000e+00   3.362e+03   4.116e+03   7.795e-02
     9  -0.00000e+00   3.150e+03   3.989e+03   7.690e-02
    10  -0.00000e+00   2.968e+03   3.884e+03   7.516e-02
    11  -0.00000e+00   2.811e+03   3.799e+03   7.289e-02
    21  -0.00000e+00   2.134e+03   3.447e+03   4.580e-02
    31  -0.00000e+00   2.026e+03   3.260e+03   3.427e-02
    32  -0.00000e+00   2.023e+03   3.242e+03   3.372e-02
    33  -0.00000e+00   2.020e+03   3.224e+03   3.322e-02
    34  -0.00000e+00   2.017e+03   3.208e+03   3.275e-02
    35  -0.00000e+00   2.015e+03   3.191e+03   3.234e-02
    36  -0.00000e+00   2.013e+03   3.175e+03   3.197e-02
    37  -0.00000e+00   2.011e+03   3.159e+03   3.162e-02
    38  -0.00000e+00   2.009e+03   3.144e+03   3.129e-02
    39  -0.00000e+00   2.008e+03   3.130e+03   3.096e-02
    40  -0.00000e+00   2.006e+03   3.116e+03   3.065e-02

FISTA finished

Linear events with subsampling in both directions

In [27]:
###### Input parameters
par = {'oy':0,'dy':2,    'ny':51,
       'ox':0,'dx':2,    'nx':30,
       'ot':0,'dt':0.004,'nt':100,
       'f0': 20, 'nfmax': 210}

v       = 1500
t0_m    = [0.1, 0.2, 0.3]
theta_m = [0, 30, -60]
phi_m   = [0, 50, 30]
amp_m   = [1., -2, 0.5]

# Create axis
taxis, taxis2, xaxis, yaxis = makeaxis(par)

# Create wavelet
wav = ricker(taxis[:41], f0=par['f0'])[0]

# Generate model
m, x =  linear3d(xaxis, yaxis, taxis, v, t0_m, theta_m, phi_m, amp_m, wav)

# subsampling locations
perc_subsampling=0.4
Nsub=int(np.round(par['ny']*perc_subsampling))
Nsub1=int(np.round(par['nx']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['ny']))[:Nsub])
iava1 = np.sort(np.random.permutation(np.arange(par['nx']))[:Nsub1])

# Create restriction operator
Rop0 = Restriction(par['ny']*par['nx']*par['nt'], iava, 
                   dims=(par['ny'], par['nx'], par['nt']), dir=0, dtype='float64')
Rop1 = Restriction(Nsub*par['nx']*par['nt'], iava1, 
                   dims=(Nsub, par['nx'], par['nt']), dir=1, dtype='float64')
Rop1mask = Restriction(par['ny']*par['nx']*par['nt'], iava1, 
                       dims=(par['ny'], par['nx'], par['nt']), dir=1, dtype='float64')

Rop  = Rop1*Rop0
dottest(Rop,Nsub*Nsub1*par['nt'], par['ny']*par['nx']*par['nt'], verb=True)

y = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, Nsub1, par['nt'])
xadj = xadj.reshape( par['ny'], par['nx'], par['nt'])

# Apply mask
ymask = Rop1mask.mask(Rop0.mask(x.flatten()))

Dot test passed, v^T(Opu)=3.886508 - u^T(Op^Tv)=3.886508
In [28]:
# Invert for interpolated signal
xinv, _, _ = SeismicInterpolation(y, (par['ny'], par['nx']), iava, iava1=iava1, kind='spatial', dottest=True,
                                  **dict(epsRs=[np.sqrt(0.1)], returninfo=False,
                                         damp=np.sqrt(1e-4), iter_lim=100, show=0))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[par['ny']//2].T-xinv[par['ny']//2].T, cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[:, :, par['nt']//2], cmap='gray', vmin=-2, vmax=2)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[:, :, par['nt']//2], cmap='gray', vmin=-2, vmax=2)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[:, :, par['nt']//2], cmap='gray', vmin=-2, vmax=2)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[:, :, par['nt']//2], cmap='gray', vmin=-2, vmax=2)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[:, :, par['nt']//2]-xinv[:, :, par['nt']//2], cmap='gray', vmin=-2, vmax=2)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=201.611621 - u^T(Op^Tv)=201.611621

Realistic data

In [29]:
inputfile3d = '../data/seismicinterpolation/input3d.npz'
inputdata3d = np.load(inputfile3d)

x = np.transpose(inputdata3d['data'], (2, 1, 0))[:,:, :400:2]
x = x/np.abs(x).max()

taxis, xaxis, yaxis = inputdata3d['t'][:400:2], inputdata3d['rx'][:], inputdata3d['ry'][:]
par = {}
par['ny'], par['nx'], par['nt'] = x.shape
par['dx'] = inputdata3d['rx'][1] - inputdata3d['rx'][0]
par['dy'] = inputdata3d['ry'][1] - inputdata3d['ry'][0]
par['dt'] = inputdata3d['t'][1] - inputdata3d['t'][0]

# add wavelet
wav3d = wav[::3]
wav3d_c = np.argmax(wav3d)
x = np.apply_along_axis(convolve, 2, x, wav3d, mode='full')
x = x[:, :, wav3d_c:][:, :, :par['nt']]

# gain
gain = np.transpose(np.tile((taxis**2)[:, np.newaxis, np.newaxis], (1, par['ny'], par['nx'])), (1, 2, 0))
x = x*gain

# subsampling locations
perc_subsampling=0.6
Nsub=int(np.round(par['ny']*perc_subsampling))

iava = np.sort(np.random.permutation(np.arange(par['ny']))[:Nsub])

# Create restriction operator
Rop = Restriction(par['ny']*par['nx']*par['nt'], iava, 
                  dims=(par['ny'], par['nx'], par['nt']), dir=0, dtype='float64')
dottest(Rop,Nsub*par['nx']*par['nt'], par['ny']*par['nx']*par['nt'], verb=True)

y = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, par['nx'], par['nt'])
xadj = xadj.reshape( par['ny'], par['nx'], par['nt'])

# Apply mask
ymask = Rop.mask(x.flatten())

Dot test passed, v^T(Opu)=640.111950 - u^T(Op^Tv)=640.111950

Inversion in time-space domain with spatial smoothing

In [30]:
# Invert for interpolated signal
xinv, _, _ = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='spatial', dottest=True,
                                  **dict(epsRs=[np.sqrt(0.1)], returninfo=False,
                                         damp=np.sqrt(1e-4), iter_lim=100, show=0))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[4].set_title('Error')
axs[4].axis('tight');

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[par['ny']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[par['ny']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xadj[par['ny']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Adjoint')
axs[2].axis('tight');
axs[3].imshow(xinv[par['ny']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Reconstructed data (reg. inv.)')
axs[3].axis('tight')
axs[4].imshow(x[par['ny']//2].T-xinv[par['ny']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[4].set_title('Error')
axs[4].axis('tight');

Dot test passed, v^T(Opu)=-798.969647 - u^T(Op^Tv)=-798.969647

Inversion in f-k domain with sparsity-promoting solver

In [31]:
nfft=2**8
FFTop = FFTND(dims=[par['ny'], par['nx'], par['nt']], nffts=[nfft, nfft, nfft], 
              sampling=[par['dy'], par['dx'], par['dt']])
dottest(FFTop, nfft*nfft*nfft, par['ny']*par['nx']*par['nt'], complexflag=2, verb=True)
p = FFTop*x.flatten()
p = np.reshape(p, (nfft, nfft, nfft))

pmask = FFTop*ymask.flatten()
pmask = np.reshape(pmask, (nfft, nfft, nfft))

# Inverse
xinv, pinv, _ = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='fk', nffts=(nfft, nfft, nfft), 
                                     sampling=(par['dy'], par['dx'], par['dt']), dottest=False,
                                     **dict(niter=50, eps=1e-4, alpha=1, 
                                            returninfo=True, show=True))

fig, axs = plt.subplots(2, 4, figsize=(20, 9))
axs[0][0].imshow(np.fft.fftshift(np.abs(pmask[:, 1, :nfft//2-1]),axes=0).T, 
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      FFTop.fs[2][nfft//2-1], FFTop.fs[2][0]), vmin=0, vmax=0.01)
axs[0][0].set_title('Masked FFT')
axs[0][0].axis('tight')
axs[0][0].set_ylim(150, 0)
axs[0][0].set_xlim(-0.02, 0.02)
axs[0][1].imshow(np.fft.fftshift(np.abs(p[:, 1, :nfft//2-1]),axes=0).T, 
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      FFTop.fs[2][nfft//2-1], FFTop.fs[2][0]), vmin=0, vmax=0.01)
axs[0][1].set_title('Original FFT')
axs[0][1].axis('tight')
axs[0][1].set_ylim(150, 0)
axs[0][1].set_xlim(-0.02, 0.02)
axs[0][2].imshow(np.fft.fftshift(np.abs(pinv[:, 1, :nfft//2-1]), axes=0).T,
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      FFTop.fs[2][nfft//2-1], FFTop.fs[2][0]), vmin=0, vmax=0.01)
axs[0][2].set_title('Reconstructed FFT')
axs[0][2].axis('tight')
axs[0][2].set_ylim(150, 0)
axs[0][2].set_xlim(-0.02, 0.02)
axs[1][0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][0].set_title('Data')
axs[1][0].axis('tight');
axs[1][1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][1].set_title('Masked data')
axs[1][1].axis('tight')
axs[1][2].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][2].set_title('Reconstructed data (sparse inv.)')
axs[1][2].axis('tight');
axs[1][3].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][3].set_title('Error')
axs[1][3].axis('tight');
                                  
fig, axs = plt.subplots(1, 6, figsize=(20, 9))
axs[0].imshow(np.fft.fftshift(np.fft.fftshift(np.abs(p[:, :, 50]),axes=0),axes=1).T, 
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      np.fft.fftshift(FFTop.fs[1])[0], np.fft.fftshift(FFTop.fs[1])[-1]), 
              vmin=0, vmax=0.01)
axs[0].set_title('Original FFT')
axs[0].axis('tight')
axs[0].set_ylim(-0.02, 0.02)
axs[0].set_xlim(-0.02, 0.02)
axs[1].imshow(np.fft.fftshift(np.fft.fftshift(np.abs(pinv[:, :, 50]),axes=0),axes=1).T, 
              extent=(np.fft.fftshift(FFTop.fs[0])[0], np.fft.fftshift(FFTop.fs[0])[-1], 
                      np.fft.fftshift(FFTop.fs[1])[0], np.fft.fftshift(FFTop.fs[1])[-1]), 
              vmin=0, vmax=0.01)
axs[1].set_title('Reconstructed FFT')
axs[1].axis('tight')
axs[1].set_ylim(-0.02, 0.02)
axs[1].set_xlim(-0.02, 0.02)
axs[2].imshow(x[:, :, par['nt']//2], cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Data')
axs[2].axis('tight');
axs[3].imshow(ymask[:, :, par['nt']//2], cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Masked data')
axs[3].axis('tight')
axs[4].imshow(xinv[:, :, par['nt']//2], cmap='gray', vmin=-0.005, vmax=0.005)
axs[4].set_title('Reconstructed data (sparse inv.)')
axs[4].axis('tight');
axs[5].imshow(x[:, :, par['nt']//2].T-xinv[:, :, par['nt']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[5].set_title('Error')
axs[5].axis('tight');

Dot test passed, v^T(Opu)=481.326296 - u^T(Op^Tv)=481.326296
FISTA optimization
-----------------------------------------------------------
The Operator Op has 1182600 rows and 16777216 cols
eps = 1.000000e-04	tol = 1.000000e-10	niter = 50
alpha = 1.000000e+00	thresh = 5.000000e-05
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1  -0.00000e+00   5.306e-01   9.287e-01   2.513e+00
     2  -0.00000e+00   1.479e-01   6.356e-01   7.040e-01
     3  -0.00000e+00   3.606e-02   5.657e-01   4.200e-01
     4  -0.00000e+00   1.317e-02   5.531e-01   2.166e-01
     5  -0.00000e+00   1.008e-02   5.427e-01   1.236e-01
     6  -0.00000e+00   1.021e-02   5.283e-01   1.164e-01
     7  -0.00000e+00   1.094e-02   5.135e-01   1.235e-01
     8  -0.00000e+00   1.149e-02   4.999e-01   1.266e-01
     9  -0.00000e+00   1.157e-02   4.877e-01   1.284e-01
    10  -0.00000e+00   1.130e-02   4.767e-01   1.299e-01
    11  -0.00000e+00   1.094e-02   4.666e-01   1.307e-01
    21  -0.00000e+00   9.459e-03   4.012e-01   1.152e-01
    31   0.00000e+00   8.876e-03   3.735e-01   9.385e-02
    41  -0.00000e+00   8.591e-03   3.599e-01   7.785e-02
    42  -0.00000e+00   8.571e-03   3.589e-01   7.657e-02
    43  -0.00000e+00   8.553e-03   3.580e-01   7.534e-02
    44  -0.00000e+00   8.534e-03   3.572e-01   7.411e-02
    45  -0.00000e+00   8.517e-03   3.563e-01   7.293e-02
    46  -0.00000e+00   8.500e-03   3.556e-01   7.180e-02
    47  -0.00000e+00   8.485e-03   3.548e-01   7.070e-02
    48  -0.00000e+00   8.470e-03   3.541e-01   6.966e-02
    49  -0.00000e+00   8.456e-03   3.535e-01   6.862e-02
    50  -0.00000e+00   8.442e-03   3.529e-01   6.760e-02

FISTA finished

Inversion in parabolic Radon domain with sparsity-promoting solver

In [43]:
npy=npx=61
pxmax = 8e-6
px = np.linspace(-pxmax/2, pxmax, npx)
Radop = Radon3D(taxis, yaxis, xaxis, px, px, centeredh=True, kind='parabolic', 
                onthefly=True, engine='numba')
dottest(Radop, par['ny']*par['nx']*par['nt'], npy*npx*par['nt'], verb=True);

Dot test passed, v^T(Opu)=-55935.466696 - u^T(Op^Tv)=-55935.466696
In [44]:
Rop1 = Rop*Radop

# Radon on original data
xradadj = Radop.H*x.flatten()
xradadj = xradadj.reshape(npx, npx, par['nt'])

yradadj = Radop*xradadj.flatten()
yradadj = yradadj.reshape(par['ny'], par['nx'], par['nt'])

# Adjoint
xadj = Rop1.H*y.flatten()
xadj = xadj.reshape(npy, npx, par['nt'])

fig, axs = plt.subplots(1, 4, figsize=(20, 6))
axs[0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj[:, npx//2].T, cmap='gray', interpolation='sinc')
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj[:, par['nx']//2].T, cmap='gray')
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj[:, npx//2].T, cmap='gray', interpolation='sinc')
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');
In [45]:
it = 110

fig, axs= plt.subplots(1, 4, figsize=(20, 6))
axs[0].imshow(x[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(xradadj[:, :, it], cmap='gray', interpolation='sinc')
axs[1].set_title('Adjoint Radon on data')
axs[1].axis('tight')
axs[2].imshow(yradadj[:, :, it], cmap='gray')
axs[2].set_title('Modelled data from Adjoint Radon')
axs[2].axis('tight');
axs[3].imshow(xadj[:, :, it], cmap='gray', interpolation='sinc')
axs[3].set_title('Radon on Restricted data')
axs[3].axis('tight');
In [46]:
xinv, pinv, cost = SeismicInterpolation(y, (par['ny'], par['nx']), 
                                        iava, kind='radon-parabolic', spataxis=yaxis, spat1axis=xaxis, 
                                        taxis=taxis, paxis=px, p1axis=px, centeredh=True, dottest=False,
                                        **dict(niter=30, eps=1e-3, returninfo=True, show=True))

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(pinv[:, npx//2].T, cmap='gray')
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[4].set_title('Error')
axs[4].axis('tight');

FISTA optimization
-----------------------------------------------------------
The Operator Op has 1182600 rows and 744200 cols
eps = 1.000000e-03	tol = 1.000000e-10	niter = 30
alpha = 4.730998e-08	thresh = 2.365499e-11
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   0.00000e+00   3.331e+00   3.331e+00   4.229e-05
     2   0.00000e+00   3.295e+00   3.295e+00   4.157e-05
     3   0.00000e+00   3.250e+00   3.250e+00   5.240e-05
     4   0.00000e+00   3.198e+00   3.198e+00   6.240e-05
     5   0.00000e+00   3.140e+00   3.140e+00   7.166e-05
     6   0.00000e+00   3.077e+00   3.077e+00   8.021e-05
     7   0.00000e+00   3.009e+00   3.009e+00   8.807e-05
     8   0.00000e+00   2.939e+00   2.939e+00   9.525e-05
     9   0.00000e+00   2.867e+00   2.867e+00   1.018e-04
    10   0.00000e+00   2.793e+00   2.793e+00   1.077e-04
    11   0.00000e+00   2.719e+00   2.719e+00   1.130e-04
    21   0.00000e+00   2.034e+00   2.035e+00   1.422e-04
    22   0.00000e+00   1.975e+00   1.976e+00   1.435e-04
    23   0.00000e+00   1.918e+00   1.919e+00   1.446e-04
    24   0.00000e+00   1.862e+00   1.863e+00   1.455e-04
    25   0.00000e+00   1.809e+00   1.810e+00   1.463e-04
    26   0.00000e+00   1.757e+00   1.758e+00   1.470e-04
    27   0.00000e+00   1.707e+00   1.708e+00   1.475e-04
    28   0.00000e+00   1.659e+00   1.660e+00   1.479e-04
    29   0.00000e+00   1.613e+00   1.614e+00   1.482e-04
    30   0.00000e+00   1.568e+00   1.569e+00   1.484e-04

FISTA finished
In [47]:
fig, axs = plt.subplots(1, 5, figsize=(20, 9))
axs[0].imshow(pinv[:, :, it], cmap='gray', interpolation='sinc')
axs[0].set_title('L1 model in sparse domain')
axs[0].axis('tight')
axs[1].imshow(x[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Data')
axs[1].axis('tight');
axs[2].imshow(ymask[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Masked data')
axs[2].axis('tight')
axs[3].imshow(xinv[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Reconstructed data (sparse inv.)')
axs[3].axis('tight');
axs[4].imshow(x[:, :, it]-xinv[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[4].set_title('Error')
axs[4].axis('tight');

Inversion in sliding linear Radon domain with sparsity-promoting solver

In [37]:
dy = par['dy']
dx = par['dx']
nwin=(21, 19)
nwins=(7, 5)
nover=(6, 5)
npy=41
npx=41
pxmax = 5e-4
py = px = np.linspace(-pxmax, pxmax, npx)
dimsd = (par['ny'], par['nx'], par['nt'])
dims = (nwins[0]*npy, nwins[1]*npx, dimsd[2])

# sliding window radon with overlap
Op = Radon3D(taxis, np.linspace(-dy*nwin[0]//2, dy*nwin[0]//2, nwin[0]), 
             np.linspace(-dx*nwin[1]//2, dx*nwin[1]//2, nwin[1]), 
             px, px, centeredh=True,
             kind='linear', engine='numba')
Slidop = Sliding3D(Op, dims, dimsd, nwin, nover, (npy, npx), tapertype='cosine', design=True)
dottest(Slidop, np.prod(dimsd), np.prod(dims), verb=True)

# adjoint
Rop1 = Rop*Slidop

padj_fromx = Slidop.H*x.flatten()
padj_fromx = padj_fromx.reshape(nwins[0], nwins[1], npy, npx, par['nt'])

padj = Rop1.H*y.flatten()
padj = padj.reshape(nwins[0], nwins[1], npy, npx, par['nt'])

WARNING: (7,5) windows required...
WARNING: model wins - start0:[  0  41  82 123 164 205 246], end0:[ 41  82 123 164 205 246 287], start1:[  0  41  82 123 164], end1:[ 41  82 123 164 205]
WARNING: data wins - start0:[ 0 15 30 45 60 75 90], end0:[ 21  36  51  66  81  96 111], start1:[ 0 14 28 42 56], end1:[19 33 47 61 75]

Dot test passed, v^T(Opu)=-13399.181031 - u^T(Op^Tv)=-13399.181031
In [38]:
xinv, pinv, cost = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='sliding', 
                                          spataxis=yaxis, spat1axis=xaxis, 
                                          taxis=taxis, paxis=py, p1axis=px, nwins=nwins, nwin=nwin, 
                                          nover=nover, dottest=False,
                                          **dict(niter=40, alpha=1.1e-5,
                                                 eps=1e-1, returninfo=True, show=True))

fig, axs = plt.subplots(2, 4, figsize=(25, 12))
axs[0][0].imshow(padj_fromx[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-1, vmax=1)
axs[0][0].set_title('Adjoint Radon on data')
axs[0][0].axis('tight')
axs[0][1].imshow(padj[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-1, vmax=1)
axs[0][1].set_title('Adjoint Radon on subsampled data')
axs[0][1].axis('tight')
axs[0][2].imshow(pinv[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-0.001, vmax=0.001)
axs[0][2].set_title('Inverse Radon on subsampled data')
axs[0][2].axis('tight')
axs[1][0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][0].set_title('Data')
axs[1][0].axis('tight');
axs[1][1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][1].set_title('Masked data')
axs[1][1].axis('tight');
axs[1][2].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][2].set_title('Reconstructed data (sparse inv.)')
axs[1][2].axis('tight')
axs[1][3].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][3].set_title('Error')
axs[1][3].axis('tight');

FISTA optimization
-----------------------------------------------------------
The Operator Op has 1182600 rows and 11767000 cols
eps = 1.000000e-01	tol = 1.000000e-10	niter = 40
alpha = 1.100000e-05	thresh = 5.500000e-07
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1  -0.00000e+00   3.299e+00   3.323e+00   8.032e-04
     2  -0.00000e+00   3.233e+00   3.279e+00   7.820e-04
     3  -0.00000e+00   3.153e+00   3.226e+00   9.771e-04
     4  -0.00000e+00   3.060e+00   3.167e+00   1.152e-03
     5  -0.00000e+00   2.959e+00   3.102e+00   1.310e-03
     6  -0.00000e+00   2.851e+00   3.034e+00   1.450e-03
     7  -0.00000e+00   2.739e+00   2.964e+00   1.573e-03
     8  -0.00000e+00   2.624e+00   2.895e+00   1.680e-03
     9  -0.00000e+00   2.509e+00   2.826e+00   1.772e-03
    10  -0.00000e+00   2.394e+00   2.759e+00   1.849e-03
    11  -0.00000e+00   2.282e+00   2.695e+00   1.912e-03
    21  -0.00000e+00   1.403e+00   2.262e+00   1.999e-03
    31  -0.00000e+00   9.754e-01   2.109e+00   1.651e-03
    32  -0.00000e+00   9.495e-01   2.100e+00   1.614e-03
    33  -0.00000e+00   9.257e-01   2.092e+00   1.577e-03
    34  -0.00000e+00   9.039e-01   2.084e+00   1.542e-03
    35  -0.00000e+00   8.839e-01   2.077e+00   1.507e-03
    36  -0.00000e+00   8.656e-01   2.069e+00   1.474e-03
    37  -0.00000e+00   8.488e-01   2.062e+00   1.441e-03
    38  -0.00000e+00   8.334e-01   2.055e+00   1.410e-03
    39  -0.00000e+00   8.193e-01   2.049e+00   1.380e-03
    40  -0.00000e+00   8.063e-01   2.042e+00   1.352e-03

FISTA finished
In [39]:
fig, axs = plt.subplots(1, 4, figsize=(20, 9))
axs[0].imshow(x[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight');
axs[1].imshow(ymask[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xinv[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Reconstructed data (sparse inv.)')
axs[2].axis('tight');
axs[3].imshow(x[:, :, par['nt']//4]-xinv[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Error')
axs[3].axis('tight');

Finally we repeat the same exercise but with an regular coarse sampling along y axis

In [76]:
# subsampling locations
iava = np.arange(0, par['ny'], 5)
Nsub = len(iava)

# Create restriction operator
Rop = Restriction(par['ny']*par['nx']*par['nt'], iava, 
                  dims=(par['ny'], par['nx'], par['nt']), dir=0, dtype='float64')
dottest(Rop,Nsub*par['nx']*par['nt'], par['ny']*par['nx']*par['nt'], verb=True)

y = Rop*x.flatten()
xadj = Rop.H*y.flatten()
y = y.reshape(Nsub, par['nx'], par['nt'])
xadj = xadj.reshape( par['ny'], par['nx'], par['nt'])

# Apply mask
ymask = Rop.mask(x.flatten())

fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Masked data')
axs[1].axis('tight')


fig, axs = plt.subplots(1, 5, figsize=(20, 6))
axs[0].imshow(x[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight')
axs[1].imshow(ymask[:, :, it], cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Masked data')
axs[1].axis('tight')

Dot test passed, v^T(Opu)=-119.949205 - u^T(Op^Tv)=-119.949205
Out[76]:
(-0.5, 80.5, 120.5, -0.5)
In [89]:
dy = par['dy']
dx = par['dx']
nwin=(21, 19)
nwins=(7, 5)
nover=(5, 5)
npy=41
npx=41
pxmax = 5e-4
py = px = np.linspace(-pxmax, pxmax, npx)
dimsd = (par['ny'], par['nx'], par['nt'])
dims = (nwins[0]*npy, nwins[1]*npx, dimsd[2])

# sliding window radon with overlap
Op = Radon3D(taxis, np.linspace(-dy*nwin[0]//2, dy*nwin[0]//2, nwin[0]), 
             np.linspace(-dx*nwin[1]//2, dx*nwin[1]//2, nwin[1]), 
             px, px, centeredh=True,
             kind='linear', engine='numba')
Slidop = Sliding3D(Op, dims, dimsd, nwin, nover, (npy, npx), tapertype='cosine', design=True)
dottest(Slidop, np.prod(dimsd), np.prod(dims), verb=True)

# adjoint
Rop1 = Rop*Slidop

padj_fromx = Slidop.H*x.flatten()
padj_fromx = padj_fromx.reshape(nwins[0], nwins[1], npy, npx, par['nt'])

padj = Rop1.H*y.flatten()
padj = padj.reshape(nwins[0], nwins[1], npy, npx, par['nt'])

WARNING: (7,5) windows required...
WARNING: model wins - start0:[  0  41  82 123 164 205 246], end0:[ 41  82 123 164 205 246 287], start1:[  0  41  82 123 164], end1:[ 41  82 123 164 205]
WARNING: data wins - start0:[ 0 16 32 48 64 80 96], end0:[ 21  37  53  69  85 101 117], start1:[ 0 14 28 42 56], end1:[19 33 47 61 75]

Dot test passed, v^T(Opu)=34991.440704 - u^T(Op^Tv)=34991.440704
In [90]:
xinv, pinv, cost = SeismicInterpolation(y, (par['ny'], par['nx']), iava, kind='sliding', 
                                        spataxis=yaxis, spat1axis=xaxis, 
                                        taxis=taxis, paxis=py, p1axis=px, nwins=nwins, nwin=nwin, 
                                        nover=nover, dottest=False,
                                        **dict(niter=50, alpha=2.030948e-6,
                                             eps=5e-2, returninfo=True, show=True))

fig, axs = plt.subplots(2, 4, figsize=(25, 12))
axs[0][0].imshow(padj_fromx[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-1, vmax=1)
axs[0][0].set_title('Adjoint Radon on data')
axs[0][0].axis('tight')
axs[0][1].imshow(padj[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-1, vmax=1)
axs[0][1].set_title('Adjoint Radon on subsampled data')
axs[0][1].axis('tight')
axs[0][2].imshow(pinv[nwins[0]//2, :, :, npx//2].reshape(nwins[1]*npy, par['nt']).T, 
                 cmap='gray', vmin=-0.001, vmax=0.001)
axs[0][2].set_title('Inverse Radon on subsampled data')
axs[0][2].axis('tight')
axs[1][0].imshow(x[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][0].set_title('Data')
axs[1][0].axis('tight');
axs[1][1].imshow(ymask[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][1].set_title('Masked data')
axs[1][1].axis('tight');
axs[1][2].imshow(xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][2].set_title('Reconstructed data (sparse inv.)')
axs[1][2].axis('tight')
axs[1][3].imshow(x[:, par['nx']//2].T-xinv[:, par['nx']//2].T, cmap='gray', vmin=-0.005, vmax=0.005)
axs[1][3].set_title('Error')
axs[1][3].axis('tight');

FISTA optimization
-----------------------------------------------------------
The Operator Op has 405000 rows and 11767000 cols
eps = 5.000000e-02	tol = 1.000000e-10	niter = 50
alpha = 2.030948e-06	thresh = 5.077370e-08
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1  -0.00000e+00   1.168e+00   1.169e+00   7.138e-05
     2  -0.00000e+00   1.164e+00   1.167e+00   7.118e-05
     3  -0.00000e+00   1.160e+00   1.165e+00   9.097e-05
     4  -0.00000e+00   1.155e+00   1.162e+00   1.101e-04
     5  -0.00000e+00   1.150e+00   1.159e+00   1.287e-04
     6  -0.00000e+00   1.143e+00   1.155e+00   1.469e-04
     7  -0.00000e+00   1.136e+00   1.151e+00   1.646e-04
     8  -0.00000e+00   1.128e+00   1.147e+00   1.820e-04
     9  -0.00000e+00   1.119e+00   1.142e+00   1.990e-04
    10  -0.00000e+00   1.110e+00   1.137e+00   2.157e-04
    11  -0.00000e+00   1.100e+00   1.132e+00   2.319e-04
    21  -0.00000e+00   9.731e-01   1.067e+00   3.706e-04
    31  -0.00000e+00   8.219e-01   9.964e-01   4.611e-04
    41  -0.00000e+00   6.760e-01   9.400e-01   5.025e-04
    42  -0.00000e+00   6.626e-01   9.355e-01   5.042e-04
    43  -0.00000e+00   6.493e-01   9.312e-01   5.055e-04
    44  -0.00000e+00   6.364e-01   9.271e-01   5.064e-04
    45  -0.00000e+00   6.237e-01   9.233e-01   5.069e-04
    46  -0.00000e+00   6.113e-01   9.197e-01   5.070e-04
    47  -0.00000e+00   5.992e-01   9.163e-01   5.067e-04
    48  -0.00000e+00   5.873e-01   9.131e-01   5.061e-04
    49  -0.00000e+00   5.758e-01   9.102e-01   5.052e-04
    50  -0.00000e+00   5.646e-01   9.074e-01   5.040e-04

FISTA finished
In [91]:
fig, axs = plt.subplots(1, 4, figsize=(20, 9))
axs[0].imshow(x[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[0].set_title('Data')
axs[0].axis('tight');
axs[1].imshow(ymask[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[1].set_title('Masked data')
axs[1].axis('tight')
axs[2].imshow(xinv[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[2].set_title('Reconstructed data (sparse inv.)')
axs[2].axis('tight');
axs[3].imshow(x[:, :, par['nt']//4]-xinv[:, :, par['nt']//4], cmap='gray', vmin=-0.005, vmax=0.005)
axs[3].set_title('Error')
axs[3].axis('tight');
In [ ]: